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Abstract. In this paper we apply the innovative Laplace transformation 
method introduced by Sheen, Sloan, and Thomee (IMA J. Numer. Anal., 2003) 
to solve the Black-Scholcs equation. The algorithm is of arbitrary high con- 
vergence rate and naturally parallelizable. It is shown that the method is very 
efficient for calculating various option prices. Existence and uniqueness prop- 
erties of the Laplace transformed Black-Scholes equation are analyzed. Also 
a transparent boundary condition associated with the Laplace transformation 
method is proposed. Several numerical results for various options under vari- 
ous situations confirm the efficiency, convergence and parallelization property 
of the proposed scheme. 
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1. Introduction 

As stock markets have become more sophisticated, so have their products. The 
simple buy/sell trades of the early markets have been replaced by more complex 
financial options and derivatives. These contracts can give investors various oppor- 
tunities to tailor their dealings to their investment needs. 

One of the main concerns about financial options is what the exact values of 
options are. For the simplest model in the case of constant coefficients, an exact 
pricing formula was derived by Black and Scholes, known as the Black-Scholcs 
formula. However, in the general case of time and space dependent coefficients the 
exact pricing formula are not yet established, and thus numerical solutions have 
been used. 

In order to describe an option price, let x, K, t and T denote the underlying 
asset price, the strike price, the time to maturity, and the expiry date of an option, 
respectively. As usual, a and r represent the volatility of the underlying asset and 
the risk-free interest rate of the market, respectively. In this paper, we assume 
that a and r depend on x only. Then a European option price u{x,t) satisfies the 
Black-Scholes equation: 

du loo d 2 u du . . „. 

(L) ~dt~2 d^~ rx dx~ + ru= ' vM)e(o,oc)x(o,T], 
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where an initial condition Uq(x) = u(x,0) is given by the initial contract of an 
option. The basket option based on n assets x = [x\, . . . , x n ) satisfies 



where <Zjj = Y^k=i ^ik^jk, with cry representing the corelation between the assets 
Xi and Xj. 

Several numerical methods have been used for solving the Black-Scholes equation, 
for example in J53J [25] and [TO] and the references therein, one can find popular 
numerical schemes for option pricing. Usually the time marching methods such as 
forward Euler, backward Euler and Crank-Nicolson schemes are used with a suitable 
spatial discretization scheme. In spite of the popularity of these time marching 
methods, a critical drawback of these schemes is that they usually require as many 
time steps as spatial meshes to balance the errors arising from discretization. In 
particular, for the estimation of basket options of reasonable size, the usual time 
marching schemes seem to be too slow in practice since the cost of solving an 
elliptic system to advance to a next time step is usually expensive. It is thus highly 
desirable to solve as small a number of elliptic solution steps as possible as well as 
to apply a very fast elliptic solver. 

In this paper, we will focus on minimizing the number of elliptic solution steps by 
proposing the Laplace transformation method for the Black-Scholes equation, which 
is also naturally parallclizable. It will be shown that our method can dramatically 
reduce the computing time compared to the time marching schemes. Suitable 
contours should be chosen in order to have very fast convergence, and for this, we 
will estimate the resolvent of the Black-Scholes equation. Also, an exact transparent 
boundary condition will be given at which the infinite spatial domain is truncated. 

There have been some related works in which the Laplace transformation method 
has been used, for instance in [3 [l8j [25] • However, in these earlier papers the 
Laplace transformation method has been used to obtain the analytic solution of 
various options rather than to develop an efficient numerical scheme. In particular, 
in [18] the partial Laplace transformation is applied for American option pricing, 
and in [51 [23] the Mellin transformation which is similar to the Laplace transfor- 
mation is used to evaluate the analytic solution of an option. Related with Laplace 
transformation methods there are other approaches based on the so-called 7i-matrix 
approach; for instance, see [8j [9], and so on. Also, high-dimensional parabolic prob- 
lems can be solved using sparse grids [HI UH HH [27] . Application of our Laplace 
transformation method using sparse grids to option pricing will also be interesting. 
Other approaches in the fast time-stepping methods can be found in [531 HOI HH] ■ 

In the following section, we will briefly describe the Laplace transformation 
method proposed by Sheen, Sloan, and Thomee in [30] with its numerical pro- 
cedure and convergence. Then in $3] we will examine the properties of the Laplace 
transformed Black-Scholes equation including the solvability of the transformed 
equation, transparent boundary condition and the resolvent. Finally in 331 we will 
present several numerical results for various options and various situations with the 
parallelization property of the proposed scheme. 



(1.2) 




(x,t)e(0,oo) n x(0,T], 
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2. The Laplace transformation and its inversion 

We begin with the abstract setting of a parabolic type equation so that the 
proposed scheme can be applicable to various problems. Consider 

OVL 

(2.1) —+Au = f, 1 6(0.71; u(Q)=u , 

where uo is a given initial function and A a spatial elliptic operator with its eigen- 
values being located in the right half plane. (We added the source term f(x,t), 
which is not present in or in (jl.2p . in order to describe our method in more 
general setting.) For each z in the complex plane, recall that the standard Laplace 
transform in time of a function u(-,t) is given by 



u{-,z) :=C[u](z) = / u(-,t)e~ zt dt. 
Jo 

Then the Laplace transformation of (|2.ip is thus given in the form 

(2.2) zu + Au = u + f(-,z), zeT, 
from which the solution u(z) — u(-, z) is formally given by 

(2.3) u(;z) = (zl + ArVoO) + /(•,*)), 

for each z. We suppose that the real parts of singular points of f(z) are less than 
some positive number. 

The Laplace inversion formula ( 2|) is given by 

(2.4) u(-,t) = -L / u(-,z)e zt dz, 

where the integral contour T is a straight line parallel to the imaginary axis ex- 
pressed as 

(2.5) T :— {z €z <C : z(lj) = a + iw where u) € K increases from — ooto + oo}. 

The constant a £ K in the contour is called the Laplace convergence abscissa, and 
the value of a is required to be greater than the real part of any singularity of u(z). 
Inserting the explicit form of z S T given by (|2.5[) into Equation (|2.4p . one has 

e at f°° 

(2.6) u(x,t) — / [Re{w(a;, a + icu)} cos(wt) — Im{u(x, a + iuj)} sin(wt)] dto. 



Denoting by the summation with its first and the last summands being halved, 
an application of the composite trapezoidal rule to this integral leads to the direct 
method 



e at t N-l 

u(x, t) w — — > 
v. j j T k= i 



A;7ri knt kni knt. 

Re{u(x, a + — )} cos( — ) - lm{u(x, a + — )} sm(— ) 



for some sufficiently large N with the length of two mesh points ir/T. Although 
this scheme can be easily implemented, its convergence rate is slow due to the 
truncation and discretization errors. In order to approximate the integration (|2.6p 
fast and accurately, there have been numerous modifications, such as [3l El 
[23 [T71 [SSI HDl H3 El [HI HH 133 [SSI [35] and the references therein. In this 
paper, we will use the deformation of the contour introduced in [3D], which gives 
an arbitrary high-order convergence rate with a hyperbolic type deformation. 
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2.1. Deformation of contour. For a concrete mathematical analysis, we assume 
that the spectrum a(A) of A lies in a sector £5 such that 

o(A) C S 5 = {z e C : | argz| < 5, z £ 0, S £ (0, |)}, 

and the resolvent (zl + A)~ x of A satisfies 

M 

\\(zi + A)- 1 !) < forze^US, 

where B is a small circle at the origin. 

The first restriction is required to avoid the singular points of the integrand in 
(|2.4p . Since the problem (|2.i[) has a solution of the form 

(2.7) u(t) (= U (., t)) = jf (zJ + A)- 1 (uo + /(z))e zt efe, 

i/ie integral contour has to be kept away from the spectrum of —A and the singular 
points of f(z), when we deform the contour. In particular, since all eigenvalues of 
—A and the singularities of f(z) have real parts bounded by a positive number, 
this restriction is natural. 

Observe that if z £ T has negative real parts as \z\ becomes large, the dis- 
cretization error in numerically evaluating the integrand in (|2.7[) will be reduced 
for positive t; thus it will be desirable to deform the contour to the left half plane 
as long as all the singularities are to the left of it. Based on this, Sheen et al. [30] 
proposed the smooth contour of hyperbola type as follows: 

f = {zeC: z(u>) — C(w) + isu/, w£K, u> increasing}, 



where C(u>) = 7 — Vu 2 + v 2 . In this case, since the contour cuts the real line at 
7 — v, 7 and v must be selected such that 7 — v is larger than the negative of the 
smallest eigenvalue of A and the real parts of singularities of f{z). Also s should 
be chosen such that all the singularities of u(-, z) be to the left of the contour Y. 

Using the above deformed contour, the inversion formula can be written as an 
infinite integral with respect to a real variable, 

u (., t) = — / «(-,CM + *»w)(C» + wJeKW+^dw. 

27T1 J-oc 

The infinite range of the above integration can be changed into to a finite region 
by the change of variables of the form 

y(u) = tanh (^) and u{y) = - tantT 1 ^) = - log 

2 r t 1 — y 

for some r > and y £ (—1,1). The above change of variables reduces from an 
integral on an infinite interval to one on a finite interval as follows: 

(2.8) 

u(;t) = J u(; C(w(y)) + i S w(v))(C'(w(tf)) + is)e^^+ is ^ t u,'(y)dy. 

2.2. Semi-discrete approximation. The last integral formula ()2.8j) in the pre- 
vious section can be discretized in time using a quadrature rule. Explicitly the 
semi-discrete approximation of u(-,t) is given by 

(2.9) usA;t) = ±± E -(■^^y^ 

j=-N+l y 
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where 

j 

Zj = z(b)j), Uj=uj(yj) and ?/,• = — , for - N < j < N. 

It is proved in [30] that the quadrature scheme (|2.9|) is of arbitrary high-order 
spectral convergence rate if in particular the source term / has high-order regularity, 
stated as follows: 

Theorem 2.1 (Sheen-Sloan-Thomee). Let u(t) be the solution of H2.1\) and let 

UN,r{t) be its approximation defined by H2.9\) . Assume that f(z) is analytic to the 
right of the contour T and continuous onto T, with pi\z) bounded on T for j < r 
and r an integer > 1, Then, for t > rr 
(2.10) 

\\U N , T (t)-u(t)\\ < ^(l + r + l) e ^(l + log + ^)(|| Wo ||+maxsup||/( fe )(z)||). 

Three important remarks should be stressed. 

Remark 2.2. The implication of the above theorem without source term f as in 
our option pricing is that the scheme is of order 0(-^ F ) with an arbitrarily large 

r > since f = is certainly analytic and f^ r '(z) is bounded on T for positive 
integer r. This implies that the discretization errors in the time direction using the 
Laplace transformation method will be negligible compared to those caused from the 
spatial discretization part in solving the Black-Scholes equation. 

Remark 2.3. In the summand (12. 9|) . an important observation is that 



dz ^duj 

w(^)^K)-^(yj)> j = o, • • • ,n, 

are independent oft. Therefore, we only have to approximate u(-,Zj) only once 
by solving the complex-valued elliptic problem (|2.2[) for a set of Zj,j — 0, 1, • • • , N. 
Then, if we need the option pricing at a different time t, the same set of spatial 
solutions u(-,Zj),j = 0, 1, • • • , N, can be used in the evaluation of the summation 
(|2.9|) with the only change in e Zjt , for the needed time t. 



Remark 2.4. Notice that each elliptic problem (|2.2p for a Zj from the set of 
Zj,j = 0, 1, • • • , N, is independent of other elliptic problems for the remaining Zj 's. 
This will minimize communication times in solving the elliptic problems (|2.2[) in 
parallel by assigning each processor to solve an independent elliptic problem without 
communicating with other processors during solving its assigned problem. 

3. Laplace transformation method for the Black-Scholes equation 

In this section, we will apply the Laplace transformation method to the Black- 
Scholes equation depending on one stock asset. A basket option depending on 
several assets can be extended from the following numerical scheme and analyzed 
in a similar way. Taking Laplace transforms of (jl.ip . we have 

1 (9^ u du 

(3.1) zu — -o- 2 x 2 ——, - — rx— — \-ru = u Q , (i,z)gl + xr. 

2 ox z ox 

In what follows, we will examine the solvability of the above equation and the 
resolvent of the Black-Scholes equation. 



6 



H. LEE AND D. SHEEN 



3.1. The weak formulation of the Laplace transformed equation. For a 

concrete mathematical analysis, we restrict our attention to a European put op- 
tion. Since the boundary condition of a put option vanishes at infinity, the partial 
differential equation can be reformulated as a weak problem in a weighted Sobolev 
space. Let L 2 (R+) be the space of square integrable complex- valued functions on 
R + which is endowed with the inner-product (v, w) = J R+ v(x)w{x) dx and the 

norm HvU^arR.) = sj (v, v). Then following pQ, the weighted Sobolev spaces are 
defined: 

Definition 3.1. Let V be the weighted Sobolev space defined by 

dv 

V = {v G L 2 (R+) : X — £ L 2 (R+)}, 
equipped with the the semi-norm and the norm 

2 



Mv 



dv 
dx 



dx 



dv 
dx 



dx 



Similarly, let Z be the weighted Sobolev space defined by 

dv 



Z = {ve l/°°(J 



X 7T e L ° 
dx 



equipped with the the semi-norm and the norm 
dv 



\v\ z = ess. sup xeI 



dx 



v\\ z = max ess. sup^™ \v\ , ess. sup^j 



dv 
dx 



Since the boundary value vanishes at infinity, we have the following Poincare- 
type inequality, which is an extension of the real- valued version, Lemma 2.7 given 

in[I]: 



Lemma 3.2. The following bound holds: 
(3-2) IMU'CRh.) < 2|w|v 



Wv G V. 



Proof. Let v G V be arbitrary. Then, by integration by parts, the following relation 
holds: 



vv dx 



xvv x dx 



xvv x dx. 



Thus we obtain 



\v\ 2 dx<2H \v\ 2 d x yn \xv x \ 2 dx^j 



This completes the proof. 



□ 



From now on, assume that the initial data uq G V, where V is the dual space 
of V. Denote by V the topological dual space of V with the norm defined by 



|u||v 



sup , 
wev\{o} iFllv 



where (•, •) is the duality pairing of V' and V. 

Then, multiplying (|3.1[) by a test function v G V and integrating on R + , one 
obtains the weak problem of (|3.1|) as follows: For each z G T, find u(z) G V such 
that 



(3-3) 



A z (u,v) = (u ,v) Vv g V, 
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where the bilinear form A z (-, ■) : V x V — > C is defined by 
(3.4) A z (u,v) = z(u,v) + B(u, v) Vu, »eV. 

where 



n/ \ 1 / 2/ \ 2 ^ U ^ W , If 2 / \ / xCcr\ o»u_ , 

B(u,v) — - a(x)x——dx + / — r(a;j + a (x) + xa[x) — )x— v ax 
2 Jr + 02: ox J R+ V ox J ox 

+ / r(x)uvdx, 



The bilinear form v4 2 (-, •), of course, depends on z G T, and so does the solution m. 

Assumption 3.3. Assume that a E Z and r 6 L°°(M + ). Moreover, assume that 
there exists a positive constant g_ such that for all x £ R+ such that 



< g_ < a(x). 



Set 



_ | (IMU~(R+) - v 2 ) 2 /(v) 2 i if is a constant, 

\(IMU=o (R+) + 2\\cr\\ 2 z ) 2 / (a) 2 , otherwise. 

We now have the following two lemmas for the continuity and coercivity of A z (-,-) : 
VxV^C. 

Lemma 3.4. Under Assumption VS. 'A the bilinear form A z (-,-) : V x V -> C is 
continuous. 



Proof. Let u,v <E V. Then, 



r 2{x)x2 TxTx dx 



J r(x) + a 2 (x) + xa(x)*^jx^-vdx < a^/ji |«|v|M|z,=(R + ) 

< |u| v |v|v) 



+ r(x) 



' dx 



< 



(\z\ + ||r|| L =o:( R+) )|ti|y|u|y, 



where Lemma 13.21 is applied in the bound of the second inequality. Therefore the 
bilinear form A z (-, •) : V x V — > C is continuous. □ 

Lemma 3.5. Under Assumvtion \3.3[ there is a non-negative constant C\, which 
is independent of u and z, such that for all u £ V 

9 

Re{A z (u,u)} > =-|u|» - (\z\ + Ci)M|£ 2(R+) . 
Proof. Under Assumption 13. 3[ the following result is known in [I], 



(3.5) 



Rc{B(u,u)}> ^\ U \ 2 V - v\\u\\ 2 L2 
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Let u £ V be arbitrary. Then, 

2" dxdx ~~ " 2 



1 2/ \ -rdudu er 2 2 



Re{ / f — r(x) + cr 2 (x) + xcrfa;)^-)^—— uefe} < g\fii 
Jr + V &E/ 



2 



Re{ / (z + r(s)jMMda:} < (|z| + ||r|| L oc (R+) )||u||| 2(M+) , 



where Young's inequality is used in the bound of the second inequality and \x 
depends on a. A combination of these inequalities completes the lemma. □ 

The compactness of embedding L 2 (M + ) c -> V, Lemma l3Tl and Lemma [3751 imply 
that there is a unique solution in the case of European put options. We summarize 
the above results in the following theorem. 

Theorem 3.6. Suppose uq £ V'. Then, under Assumvtion \3.3\ Problem (13. 3[) has 
a unique solution u £ V. 

3.2. Resolvent of the Black-Scholes equation. In Sj2] the resolvent of a spatial 
operator is assumed to be bounded in a given sector. This assumption for the Black- 
Scholes equation will be verified in this subsection. 

Denote by R(z, —B) = (zl + B)^ 1 the resolvent of —B, so that for each / £ V', 
v = R(z, —B)f is the solution of 

(3.6) B(v,(f>)+z(v^) = (f^}, V0 £ V. 

Then we have the following lemma, which is an extension of Lemma 2.1 in |4J. 

Lemma 3.7. Under A ssumvtion 1 3. 3\ for any 9 £ (^tt,tt) there are C = C(9) > 
and k — n{9, r, a) > 0, independent of z and f , such that 

\\R(z,-B)f\\ L2{R+) < _£-_||/|| iS(R+) , forz£^ K , e , f £ L 2 (R+) 

\Z — K\ 

where = {z £ C : | arg(z — k)\ < 9}. Explicitly, the coefficients are given by 
C = (1 + |<5)(1 + S 2 ) and k= (l + ^ fi, where 5 = tan §. 

Proof. For z £ we write 

z - « = (£ + iv? =£, 2 -V 2 + 2i$7? with £ + irj £ Y, 0fi/2 , £, r) £ R, 

for any k > 0. Setting (5 = tan|, we see that 6 > 1 and |tj| < <5£. and thus the 
following inequality holds: 

e<\z-*\=e+v 2 <(i+s 2 )e- 

Set 

F = B(v,v) + z\\v\\ 2 L 2 (R+) . 
Taking the real part of F, we obtain 

ReB(v, v) + (n + e- V 2 )\Mh(R +) = ReF. 
By the inequality (|3.5p we have 

2 

(3.7) ^n 2 , + (k + e - v 2 - ri\Mh { m + ) < \n 
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By taking the imaginary part of F, we have 

Imfl(v,«)+2£7|M|£ a(R+) =lmF, 

and since lmB(v,v) = Im J ffi+ ( — r(x) + a 2 (x) + xa(x)^^jx^v dx, 

%M M\h(R + ) < \F\+<LVn\v\v\\v\\ L *(m + )- 
Multiplying by |<5 = ^tan(^0) the last estimate, we have 

V 2 \M 2 L 2 {R+) < 6\r}\ \\v\\ 2 L 2 iR+) < ^S\F\ + ljSayfJl\v\ v \\v\\ L 2 {&+) . 
Adding this to (|3.7[) . we obtain 

£\v$ + (n + e- riM^) < (i + + f Mv + ^Iklli^). 

With the choice of 

we have the following inequality 

yHv + C 2 IHIi 2(K+) <(i + ^)m 



If / e L (R + ), we take 4> = v in (|3.6p . then we have 



l| 2 ( R+) < (i + 5*) I / fvdxl < (l + ^)||/IU 2(R+ )lkllL2 (R+) , 

«/ ]R_|_ 



IV 

and therefore 

||#(*, -B)f\\ L 2(R +) < ^ 2 2 ||/||£2 (R+) < \ Z - K \ -y\\L 2 ^+)- 

This completes the proof. □ 



l + l*„„, ^ (1 + 1*)(1 + * 2 ), 

e 



From this lemma one can determine the location of a integration contour. In 
particular, if one sets the asymptotic slope of a hyperbola as s, the contour has to 
cut the real line at a point which is larger than 

tan 2 (i arctanfs)) \ 
1+ 2 2 — J »• 

In the special cast that r(x) — r and <j{x) — a are constants, k can be given by 

/ tan 2 (i arctan(s)) 
(3.8) k= 1 + — 



3.3. The transparent boundary condition. As one can see in (11.11) or (| 1 . 2[) . 

the space domain of the underlying asset of an option is an unbounded set. To 
apply a numerical scheme, one usually truncates the infinite domain into a finite 
one, and then imposes a suitable boundary condition on the boundary. Let L be a 
sufficiently large asset price. One then has the following version of the Black-Scholcs 
equation truncated at x = L. 

(3.9) ^-IaV| 2 |-rxg + ™ = 0, (x, t) G (0, L) x (0, T], 

(3.10) u(x,t) = g(x,t), (x,t) G 9(0, Z) x (0,T], 

(3.11) u{x,0) = u Q , x€[0,L]. 
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In many cases, the boundary condition on the artificial boundary x = L is im- 
posed by extending a given payoff function. For example, European put options 
assume u(L,t) = and European call options assume = 1. In [12], the 

errors caused by Dirichlet boundary conditions on the artificial boundary are esti- 
mated and thus one can determine a suitable truncation asset price for the artificial 
boundary to meet a given error tolerance. 

Instead of such artificial boundary conditions, a transparent boundary condition 
is introduced in [1] with which one can evaluate the solution in the truncated domain 
without any truncation error. However, the boundary condition in [1] is an integro- 
differential one, which needs some suitable numerical schemes to approximate it 
that will produce other possibly significant errors. We will analyze the transparent 
boundary condition in more detail and then depart from such an integro-differential 
type, by implementing the boundary condition in the Laplace transformed setting 
instead of the usual space-time setting. Our transparent boundary condition is 
motivated by the following proposition. 



Proposition 3.8. Assume that the coefficients a and r in (|3.1j) are constants and 
that L > is sufficiently large so that supp(uo) C [0, L). Then among the solutions 
u(x,z) satisfying (|3.1|) there is a component, u + , satisfying the following: 



(3.12) ^(^z^-Lj-^-Ia 2 )- \J(r-^) 2 + 2^(r + z)ju + ( Xl z) 

Vx G (L, oo), where Rc{ tfz\ > for nonzero z G C. 

Proof. Take the change of variables, y = logx, to (|3.ip . Denoting by v its solution, 
owing to supp(ito) € [0, L), one observes that v satisfies the right exterior problem 

(3.13) z£-I CT 2 0-(r-ia 2 )^+r£=O, (y, z) G (L, oo) x V. 

Among the two linearly independent solutions, we take the component which van- 
ishes at infinity, which is given as follows: 

v+(y,z) = cxp ( ,, ., v ., 

V| o~ a z V 2 

Restoring the change of variable, x = e y , and denoting by u+(x) = v+(y), one gets 

u+(x, z) = x { ) . 

Thus, by differentiating with respect to x, one arrives at 




^±(x, z) = -L J _( r _ a y 2 ) - \J (r - + 2a\r + z) j u + (x, z). 

Thus, u + satisfies the equation (|3.12[) . which completes the proof. □ 



Due to Proposition 13.81 by choosing L > sufficiently large so that supp(wo) G 
[0, L), we propose the following transparent boundary condition at x = L : 

(3-14) 

g(L, z) = ± j -(r - \a 2 ) - + sj{r~ + 2a*(r + z) \ u(L, z) Vz G T. 
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Remark 3.9. By the Laplace inversion of (|3.14[) , the transparent boundary condi- 
tion in the space-time domain is given by 

/ _2 /2) 2 ______ 

where r\ = ^ / - — |- r. /n i/ie derivation of (13 . 1 5[) . ifte following equalities are 
used: 



z£{i}£{u(t)e"} 



— 77). 

In solving the partial integro- differential equation (|3.12p tuit/i (|3.15[) using a Crank- 
Nicolson type of time-marching algorithm, one usually needs an expensive algo- 
rithm in computing time and memory. We will compare our Laplace transformation 
method with the Crank- Nicols on method in and conclude superiority in using 
our method. 

4. Numerical results 

We applied the Laplace transformation method for time discretization while the 
standard piecewise linear (Pi ) finite element method for the space discretization is 
used. Using an analytic solution for the first two examples, we can compare the 
convergence rate of the proposed scheme. In Example 14.21 we examine the effects 
of the Dirichlct boundary condition and the transparent boundary condition (|3.14[) 
in the calculation of option prices. 

In calculating the numerical values of the analytical solution, the error function 
er/(x) is evaluated by using the algorithm on page 213 of Numerical Recipes in 
Fortran [____] which has 16-digit precision. The reduction rate and speedup are 
defined by 

, ,. , II^Ax — M eX act||L 2 (0,L) 

reduction rate = log 2 — -, 

\\U A_ — Wcxact||L2( 0i L) 

where ua-x denotes the numerical solution with the spatial mesh size Ax, and 

time consumption 

speed up = 



time consumption using 1-CPU 



Example 4.1 (European put option with constant coefficients). We consider an 
European put option with coefficients r = 0.05, cr — 0.3, T = 1.0 and K = 50 and 
we truncate the domain at L = 200. 

For the numerical solutions, the boundary condition at x = in (|3.10p is given 

by 

u(x, t) = Ke~ rt , (x, t) G {0} x (0, T], 

while that at x = L 

(4.1) u(x,t)=0, (x, t) £ {L} x (0, T]. 

Although an analytic solution to this example is given by Black and Scholes, it is 
our aim to compare convergence rates for the proposed scheme and the standard 
time-marching algorithms such as Crank- Nicolson scheme. Table [_] shows conver- 
gence rate for the Crank-Nicolson scheme. As can be expected, it gives first-order 
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convergence rate. Table [3] shows that the choice of 15 z-points in the contour with 
the proposed method is enough to obtain the same level of tolerance attained using 
640 time steps with the Crank-Nicolson method. Observe that for each z-point the 
cost of solving the complex-valued elliptic problem using the proposed method is 
almost comparable to that of advancing one step forward by solving the real- valued 
elliptic problem with the time-marching algorithms. 

In the proposed scheme, we need the value of re as in Lemma 13.71 to determine 
the location of a integration contour. Since the coefficients are constants, if we 
choose the asymptotic slope of the contour as 0.4, we have re = 0.01811 by (|3.8[) . 
and therefore the contour has to cut the real line at a point greater than 0.01811. 
Under this constraint, we choose the optimal parameters which are suggested in |37j . 
and these parameters are attached in Table [3] in the case that the evaluation time 
is 1.0 for different iteration numbers. In particular, Table |3] says that 12 iterations 
are enough to balance with the space discretization of 2560 spatial meshes. 



Time steps 


Number of 
space meshes 


Mesh size 


Error in L A 


Reduction rate 


10 


10 


20 


2.928 




20 


20 


10 


0.7536 


1.958 


40 


40 


5 


0.1878 


2.004 


80 


80 


2.5 


0.4695E-01 


2.000 


160 


160 


1.25 


0.1174E-01 


2.000 


320 


320 


5/8 


0.2934E-02 


2.000 


640 


640 


5/16 


0.7337E-03 


2.000 



Table 1. Example 14. II with the Crank-Nicolson method 



Number of z 


Number of 
space meshes 


Mesh size 


Error in L l 


Reduction rate 


15 


10 


20 


2.924 




15 


20 


10 


0.7524 


1.959 


15 


40 


5 


0.1876 


2.004 


15 


80 


2.5 


0.4688E-01 


2.000 


15 


160 


1.25 


0.1172E-01 


2.000 


15 


320 


5/8 


0.2930E-02 


2.000 


15 


640 


5/16 


0.7327E-03 


2.000 



Table 2 . Example 14.11 with the Laplace transformation method 



Example 4.2 (European put option with transparent boundary condition). We 

consider a European put option with coefficients r — 0.05, a = 0.3, T — 1.0 and 
K = 50 and we truncate the domain at L = 50. 

In this example, we truncate the domain at the strike price, and then we replace 
the Dirichlet boundary condition (|4.1| with the transparent boundary condition 
given in (|3.14[) . An identical contour as in the previous example has been adopted. 
Table 2] shows that the Dirichlet boundary condition with the domain truncation 
makes a significant error, which cannot be overcome by mesh refinement. Table [5l 
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Number 
of z 


Number 
of space 
meshes 


L 2 -Error 


Reduction 
rate 


7 


V 


s 


T 


3 


2560 


0.6397E-00 




13.48 


12.42 


0.4213 


0.16500 


6 


2560 


0.1705E-01 


5.229 


26.95 


24.84 


0.4213 


0.09385 


9 


2560 


0.3434E-03 


5.634 


40.43 


37.26 


0.4213 


0.06809 


12 


2560 


0.5642E-04 


2.605 


53.90 


49.68 


0.4213 


0.05430 


15 


2560 


0.4731E-04 


0.003 


67.38 


62.09 


0.4213 


0.04556 


18 


2560 


0.4721E-04 


0.001 


80.86 


74.51 


0.4213 


0.03947 


21 


2560 


0.4717E-04 


0.000 


94.33 


86.93 


0.4213 


0.03494 



Table 3. Contour Parameters for Example 14.11 



however, gives second order convergence which is shown in Table although its 
domain is much smaller than that for Example 14.11 Indeed, comparing the same 
mesh sizes in Table [5] and Table [2j one can observe the numerical values are almost 
identical. In Figure [T] we can see the difference between the transparent boundary 
condition and the Dirichlet boundary. 



Number of z 


Number of 
space meshes 


Mesh size 


Error in L 2 


Reduction rate 


15 


10 


5 


10.35 




15 


20 


2.5 


10.40 


-0.007 


15 


40 


1.25 


10.41 


-0.002 


15 


80 


5/8 


10.42 


0.000 


15 


160 


5/16 


10.42 


0.000 


15 


320 


5/32 


10.42 


0.000 


15 


640 


5/64 


10.42 


0.000 



Table 4. Example 14.21 with the Dirichlet boundary condition at 
L = 50 



Number of z 


Number of 
space meshes 


Mesh size 


Error in L 2 


Reduction rate 


15 


10 


5 


0.1870 




15 


20 


2.5 


0.4656E-01 


2.006 


15 


40 


1.25 


0.1163E-01 


2.001 


15 


80 


5/8 


0.2907E-02 


2.000 


15 


160 


5/16 


0.7267E-03 


2.000 


15 


320 


5/32 


0.1817E-03 


1.999 


15 


640 


5/64 


0.4551E-04 


1.998 



Table 5 . Example 14.21 with the transparent boundary condition 
dSHH) at L = 50 
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Example 4.3 (Basket option with two underlying assets). We consider a European 
put basket option with two underlying assets having coefficients r — 0.05, an = 0.09, 
a 22 = 0.09, a\2 — a<x\ — —0.018, time to maturity =1.0, artificial boundary h\ = 
300,7^2 — 300 and payoff function (100 — max(xi, 2^2))+ * s given. 

For numerical computation, the boundary conditions are given by 

— (x,t) = 0, for (x,t) £ ({0}x (0,L 2 )U(0,Li) x {0}) x (0,T], 
(4.2) u(x, t) = 0, for (x, t) g ({i x } x (0, L 2 ) U (0, L x ) x {L 2 }) x (0, T\. 

To evaluate the convergence rates for the proposed scheme, we solve the same 
problem using the Crank-Nicolson scheme on a 512 x 512 space grid for the extended 
artificial domain L\ = L2 = 600 with At = 0.02. We set this as the reference 
solution and calculate the relative L 2 error for the proposed scheme. The integration 
contour is built using the parameters 7 = 35.94, v = 33.12, s = 0.4213, r = 0.07472. 
Numerical results in Table [S] show an almost second-order convergence rate. 



Number of z 


Number of 
space meshes 


Mesh size 


Relative error 
ml? 


Reduction rate 


15 


16 x 16 


75/4 


0.3662E-01 




15 


32 x 32 


75/8 


0.1047E-01 


1.806 


15 


64 x 64 


75/16 


0.2969E-02 


1.819 


15 


128 x 128 


75/32 


0.8444E-03 


1.814 



Table 6. Convergence rate in Example 14.31 on the domain 
[0, 300] x [0, 300] 



To shorten the artificial boundary, we apply the transparent boundary condition 
by assuming that the tangential derivative is negligible on the boundary. Then the 
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boundary condition 



is replaced with 



^-(L 1 ,x 2 ,z) = — - — < -(r - \ciii) - t (r - ^«u) 2 + 2an(r + z) \ u(L 1 ,x 2 ,z) 
oxi L\a\\ I 2 V 2 

on [x\, X2,z) G {L\} x (0, L 2 ) x F, and 



du 1 

— (x 1 ,L 2 ,z) = 

ax 2 L 2 a 22 



-(r - ^a 22 ) - 



(^—-022) + 2a 22 (r + z) > u(x, L 2 , z) 



on (x\,x 2l z) e (0, Li) x {£2} x T. In Table [7] we compare the results produced 
by the different boundary conditions on the lines L\ = 150 and L 2 = 150. As 
can be seen in Table [3 the transparent boundary condition is more accurate than 
the Dirichlet boundary condition. Furthermore, Table [5] and Table [7] show that 
if we apply the transparent boundary condition, it gives competitive error level 
in comparison to the Dirichlet boundary condition even though its computational 
domain is a quarter size of that with the Dirichlet boundary conditions applied. 



Num- 
ber of 

z 


Number of 
space meshes 


Mesh size 


Relative error in 
L 2 (Dirichlet) 


Relative error in 
L 2 (Transparent) 


15 


16 x 16 


75/8 


0.1998E-01 


0.1076E-01 


15 


32 x 32 


75/16 


0.1176E-01 


0.3485E-02 


15 


64 x 64 


75/32 


0.9283E-02 


0.1724E-02 



Table 7. Effect of boundary conditions in Example 14.31 on the 
domain [0, 150] x [0, 150] 



Since the elliptic equations in (|3.1| for z = zu, k = 0, 1, 2, • • • ,N, are independent 
each other, no communication is required during the computation except for the 
last summation step in the numerical Laplace inversion. Thus the Laplace transfor- 
mation method is very well fitted for parallel computation. The result in Table [8] is 
generated on 128 x 128 space grid for L\ = L 2 = 300 with a 15-number of z points 
using IBM PowerPC97 with 2.2GHz clock speed. This table, as can be expected, 
shows almost ideal speedup because of the minimization of communication time. 
Finally, we attach the plot of the basket option price at Figure 



Number of CPUs 


1 


3 


5 


15 


Time (sec) 


74.93 


25.25 


15.31 


5.671 


Speedup 


1.00 


2.97 


4.89 


13.2 



Table 8. Parallelization speedup in Example 14.31 
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